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ABSTRACT 

Recent observations support the idea that nuclear black holes grew by gas accretion while shining as luminous 
quasars at high redshift, and they establish a relation of the black hole mass with the host galaxy's spheroidal 
stellar system. We develop an analytic model to calculate the expected impact of mergers on the masses of 
black holes in massive clusters of galaxies. We use the extended Press-Schechter formalism to generate Monte 
Carlo merger histories of halos with a mass 10 A -1 Af0. We assume that the black hole mass function at z = 2 
is similar to that inferred from observations at z = (since quasar activity declines markedly at z < 2), and 
we assign black holes to the progenitor halos assuming a monotonic relation between halo mass and black 
hole mass. We follow the dynamical evolution of subhalos within larger halos, allowing for tidal stripping, the 
loss of orbital energy by dynamical friction, and random orbital perturbations in gravitational encounters with 
subhalos, and we assume that most black holes will efficiently merge after their host galaxies (represented as 
the nuclei of subhalos in our model) merge. Our analytic model reproduces numerical estimates of the subhalo 
mass function. We find that mergers can increase the mass of the most massive black holes in massive clusters 
typically by a factor ~ 2, after gas accretion has stopped. In our ten realizations of 10 15 7i -1 Mq clusters, the 
highest initial (z = 2) black hole masses are 5 - 7 x 10 9 M©, but four of the clusters contain black holes in the 
range 1 - 1 .5 x 10 1() Mq at z = 0. Satellite galaxies may host black holes whose mass is comparable to, or even 
greater than, that of the central galaxy. Thus, black hole mergers can significantly extend the very high end of 
the black hole mass function. 

Subject headings: black hole physics — cosmology: theory — dark matter — galaxies: clusters: general — 
methods: numerical — quasars: general 



1. INTRODUCTION 

The idea that quasars are powered by a ccretion of gas 
in a disk orbiting a massive black hole (ISalpeterl Il964t 
IZel'dovich & Novikovll964tlLvnden-BellfT969l) ' requires that 
enough black holes should exist at the present time to account 
for all the mass that had to be accreted over the past history 
of the universe by all observed quasars (Sol tan|[T982h . Much 
progress has been made over the last decade by discovering 
black holes in galactic nuclei, measuring their mass es and 
characterizing the population dMagorrian et al.1 [19981). Us- 
ing the surprisingly tight relation that is observed between the 
black hole mass and the velocity dispersi on of the spheroidal 
stella r population of the host galaxy (iFerrarese & Merrittl 
I2000t iGebhardt et all 120001: iTremaine et alj 120021 and refer- 
ences therein), one finds that the average mass density of the 
population of nuclear black holes is roughly equal to what is 
required for the amount of light that has been emitted by all 
Active Galactic Nuclei (AGN), if the radiative efficiency is 
close to 10%, a typical va lue expected for geometrically thin 
accre t ion disks ("see. e.g..|Barger et aLll2001t lYu & Tremaind 
l2002t IShankar et al.ll2006l) . This implies (with the caveat of 
the substantial observational uncertainties that are still present 
in the determination of the black hole mass density and the 
integrated AGN emission) that a major contribution to the 
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growth of the black hole mass density is the accretion of 
gas that is responsible for AGN, and that contributions from 
radiatively inefficient accretion (which may arise in super- 
Eddington or strongly sub-Eddington accretion flows) are at 
most comparable to the accretion from geometrically thin 
disks. 

Black hole mergers do not change the total mass density 
in the black hole population, so they do not alter this argu- 
ment, but they can alter the mass distribution of black holes. 
Observationally, we have relatively little handle on the black 
hole mass function outside the range 10 7 - 10 9 M Q that dom- 
inates the average mass density. The most luminous quasars 
of which we have evidence a t high redshift im ply black hole 
masses close to 10 10 M Q (e.g., [Fan et a"T1l20 03). assuming they 
are radiating at the Eddington luminosity. The most mas- 
sive black hole discovered so far is in M87, with a mass of 
M. ~ 3 x 10 9 M Q . It is natural to assume that the highest 
mass black holes may be lurking in some of the most massive 
elliptical galaxies, but we do not know whether the power- 
law M. - a correlation extends up to these largest black hole 
masses. Similar uncertainties plague the determination of the 
abundance of low-mass black holes. 

A number of interesting questions can be raised in rela- 
tion to theoretical expectations for the most massive black 
holes. Under the assumption that the Eddington luminosity 
cannot be exceeded (a point recently questioned by Begel- 
man, Volonteri & Rees [2006]), what is the abundance of the 
most massive black holes we should expect? Do the cold dark 
matter models for structure formation predict that the most 
massive black holes should generally reside in central cluster 
galaxies, or can they sometimes reside in a massive galaxy 
that has only recently merged in a cluster, still moving in the 
cluster outskirts? How much could the most massive black 
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holes have grown by mergers? If we search today for black 
holes in the central cluster galaxies of the most massive clus- 
ters, what is the black hole mass we should expect to measure? 

With these questions in mind, in this paper we develop a 
dynamical model for the evolution of the supermassive black 
hole population. We use this model to examine the effect of 
black hole mergers on the high-mass end of the present-day 
mass distribution. The model we develop aims to treat care- 
fully the dynamics of a halo containing a central black hole 
once it becomes an orbiting satellite within a larger halo. The 
problem of the merger of two central black holes after their 
halos are merged will not be considered here, and we simply 
obtain the maximum merger rate by assuming that two black 
holes will always merge on a short time scale once they be- 
come bound to each other in the center of two merging halos. 
After the host halo of a black hole's parent galaxy merges 
into a bigger halo, it becomes an orbiting subhalo. The sub- 
halo may either remain in orbit indefinitely (in which case 
the black hole will not merge), or it may be dragged by dy- 
namical friction to the center of the larger halo within which 
it is orbiting and be tidally disrupted. In the latter case, a 
black hole merger may ensue. We will present the result 
of a model where black holes are initially distributed among 
host halos after they have been formed by gas accretion, and 
then their mergers are followed according to the dynamical 
evolution of their subhalos. Our m odel is similar to that of 
IVolonteri. Haardt & Madaul (120031) . though we treat the dy- 
namical evolution of the subhalos in greater detail. We shall 
focus in this paper on the black holes that are present in a 
halo of 10 15 Ii~ 1 Mq at the present time, representing a massive 
galaxy cluster. 

Our implementation of a merger tree for dark matter halos is 
described in §|2] and the method to follow the dynamical evo- 
lution of merged subhalos is explained in §[3] and §0] In §|5J 
we specify the initial conditions by which we populate halos 
with black holes of different masses at an initial redshift. The 
results on the importance of mergers for black holes in clus- 
ters, and the masses of the central black holes, are presented 
in §[6] and discussed in §|7] 

2. PROPERTIES OF DARK MATTER HALOS 

We develop in this paper a dynamical model of substructure 
in dark matter halos. The first step is to construct a merger 
tree, starting from a present day halo (which for the applica- 
tion explored in this paper will represent a halo like that of 
the Coma cluster) and going backwards in time to generate 
all the merger events. This section describes how the merger 
tree is generated, how each halo in the tree is assigned a mass 
and a radius, and the choice we make for the halo density pro- 
files. In § |3] we describe the second step, where we follow 
the dynamical evolution of the halos after they have become 
satellites of a larger halo, starting at the earliest time when the 
first low-mass halos form, and moving forward to the present 
time. 6 



of a halo of mass comparable to the M ilky Way galaxy 
(e.g., iBullock. Kravtsov & We inberg 2000), although in de- 
tail there are deviations of the halo mass function from the re- 
sults of numerical simulations (e.gjEke. Cole & Frenkll 19961 : 
ISheth & Tormenlll999tlJenkins et alJ^OoUT " 

We define n(M, z) as the number density of halos of mass 
greater than M at redshift z, and <j(M,z) as the extrapolated 
linear rms fluctuation in spheres containing an average mass 
M. The halo number density is related to the fraction of the 
mass that is in halos of mass M per unit In a(M) by 



M dn(M,z) 
p m din a 



(1) 



where p m is the mean matter density of the universe. The 
Press-Schechter model assumes that 



f=\ eX P 

V 7r a \ 2<7 Z 



(2) 



where S c (z) is the critical threshold on the linear overden- 
sity for collapse at z- We adopt the approximate fitting 
formula for 6 c (z) of Eke, Co le, & Frenk (1996; see also, 
ICarroll. Press & Turned [l 9921) for a ACDM universe. The 
idea for the merger tree is based on the conditional proba- 
bility that a particle in a halo of mass M 2 at redshift z 2 was 
part of a halo of mass Af 1 at redshift z\{> Z2) dLacev & Cold 
ll993tlCole etafll2000h . 



fi 2 (M u M 2 )dM l = 
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where each subscript indicates the corresponding redshift. 
Then, the mean number of halos of mass M\ that merge 
into a halo of mass M 2 over a time step dt\ is dN jdM\ = 
(df\ 2 /dt\)(M 2 /M\)dt\. The merger tree with mass resolution 
M res is built by generating new progenitors of m ass M re5 < 
Mi < M 2 /2 with total probability P given by dCole et alj 
2000) 



P = 



m 2 /2 



M ri , 



dN 

1m 



■dM, 



(4) 



and adding a rate of smooth mass accretion, F, equal to 



F = 



dN Mi 
dM\ M 2 



dM x , 



(5) 



2.1. Merger History 

We use the extende d Press-Schechter formal i sm to gen- 
erate the merger tree dPress & Schechterl Il974t iBond et all 
1991). This formalism is very simple to use and is known 
to provide an adequate description for a merger history 

6 Throughout the paper, we interchangeably use satellite, subhalo, and 
substructure to refer to a distinct gravitationally self-bound halo in a larger 
dark matter halo. 



to account for the mass of unresolved halos with M\ < M les 
that are not discretely generated. 

2.2. Density Profile 

Density profiles of halos from nu merical simulations have 
been fit ted to a var ie ty of fo rms (|Navarro. Frenk & White! 
19971; iMoore et al l 119991; iFukushige & Makinol |200lF 
Navarro et al. I 12004 . Here we will be interested in 
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describing profiles of satellite halos with finite mass 
aft er they have merged in to a la rger halo. The profile 
of [Navarro, Fr enk & White! (1 19971) has a density slope 
(ilogp/iilogr approaching 3 at large radius, giving rise to a 
logarithmic divergence of the total mass, which demands the 
introduction of a cutoff radius. We adopt instead the simple 
form of the Jaffe sphere: 



ph(r) = 



( M h 



\4irRlJ r 2 (r+R h ) 2 



(6) 



where the Jaffe radius Rh is defined as the radius inside which 
the mass is half of the total halo mass Mh. This profile is 
steeper than that found in purely collisionless numerical sim- 
ulations in the central part, but it has the advantage of describ- 
ing a finite mass halo with a very simple form for the potential. 
In addition, black holes should in reality be in the center of a 
galaxy, and the baryonic component will steepen the profile 
near the center making it closer to that of a Jaffe sphere. 

In later sections, we will be using the one-dimensional 
velocity dispersion of the Jaffe sphere at each radius for 
isotro pic orbits, which is found from the spherical Jeans equa- 
tion ( Binn ev & Trem aine 1987) to be 



<J 2 (r)-- 



6GM h 
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GM h 




2R h 
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2.3. Radii of Dark Matter Halos 

The evolution of a halo of mass Mh after it merges with 
another object depends not only on its mass, but on its half- 
mass radius Rh as well, because in general denser objects will 
be less vulnerable to tidal disruption and will more easily sur- 
vive as orbiting satellites. It is therefore important to esti- 
mate the radius of a halo as it evolves, first by acquiring new 
mass through mergers and accretion, and then by losing mass 
through tidal disruption once it has been incorporated into a 
larger halo as a satellite. 

When a halo first appears in the merger tree at the res- 
olution mass M re s. we assign an initial radius to the halo, 
assuming the standard collapse model in which a spherical 
top-hat overdense region turns around from the Hubble ex- 
pansion at radius R t , and then virializes at a collapse ra- 
dius R co ] = R t /2. We use the fitting formula of Bryan & 
Norman (1998) for the mean density of the non-linear col- 
lapsed object, p co i = Pcrit(^)A co i, where p cl -i t is the critical den- 
sity at the redshift z of collapse, the mean matter density is 
p m {z) = n m (z)pcrit, and a flat cosmological model is assumed: 

A co i = 1 8tt 2 + 82jc - 39x 2 , x = fJ m (z) - 1 . (8) 
The collapse radius is obtained from 



4-7T 

M b =—R 3 colPcol . (9) 



For the Jaffe density profile of halos we will be assuming in 
this paper, the Jaffe radius Rh that corresponds to this collapse 
radius is obtained by equating the initial potential energy of 
a homogeneous sphere at the turnaround radius R t = 2R co \ to 
half of the potential energy of the Jaffe sphere, which gives 



5 5 

Rh= T Real = T^Rt 

o 12 



(10) 



However, as halos continue to merge, their radius should 
depend on the detailed merger history. Halos found at a red- 
shift z in which most of the mass has been recently accreted 
should have a mean density of the order of p m \, but halos that 
acquired most of their mass at an epoch much earlier than red- 
shift z and have since then accreted at a low rate should have 
a much higher density, reflecting th e density of the universe 
at the time they were formed (e.g., iNavarro. Frenk & White 
1 1997b . To simulate this process, the radii of halos in the 
merger tree are evolved as follows: At any time step in which 
two halos of masses Mm and Mh2 merge, and in addition 
an unresolved mass M acc is accreted, we consider that the 
new halo of mass Mh = Mhi +Mh2+M acc has formed from 
a configuration where the two halos and the smoothed mass 
M acc turned around from a radius R t = 2R co i, with R co i given 
by equation ( fTOb . The total energy of the final halo with a 
Jaffe sphere profile after virialization is half its potential en- 
ergy, or -GM^/(4/?h), where Rh is the Jaffe radius of the final 
halo, and we equate this to the total energy of the system at 
turnaround: 



4R h 



4R M 



GMj 2 
4fl h2 



(11) 



3 G(M a 2 cc + 2M h iM h2 + 2M hl M acc + 2M h2 M acc ) 



2R 



col 



In general, the gravitational potential energy of a ho- 
mogeneous sphere of mass Mh and radius 2R co \ is 
(3/5)GMl/(2R col ). In the last term of equation (O, we have 
included the self-gravity of mass M acc , and the interaction 
terms of mass Mhi with Mj,2, of mass Mhi with M acc , and of 
mass Mh2 and M acc . The self-gravity terms of Mhi and Mh2 are 
not included because these are already part of the total energy 
of the two merging halos that were previously virialized. In 
the case of a time step in which there is only smooth accretion, 
we simply put Mh2 = in equation (TT2l) . 

3. SUBSTRUCTURE DYNAMICS 

This section describes the dynamical evolution of a dark 
matter halo once it merges into a larger halo within which it 
orbits as a satellite. At every point in the merger tree when two 
halos merge, the smallest of the two halos becomes a satellite, 
and the large one is the host halo within which the satellite 
orbits. The mass of a halo increases at every step due to ac- 
cretion and mergers until the moment it merges with a halo 
larger than itself. After this time, the mass of the halo as a 
satellite can only decrease due to tidal stripping. 

The model we use to compute the dynamical evolu- 
tion is based on ideas sim i lar to those in the models of 
iTavlor & Babul d200ll[2003bJBullock Kravtsov & Weinberg! 
( 2000. 1200 lb. and IZentner & Bullock! (12003b. The major dif- 
ferences are the following: (1) We compute the dynami- 
cal evolution of all the halos identified in the merger tree, 
instead of following only the history of the most massive 
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halo that eventually becomes the main halo at the present 
day. (2) We continue to follow the evolution of satel- 
lites after their host halo becomes a satellite itself. (3) 
We include the random variations of orbits due to relax- 
ation of the satellites with each other. This type of analytic 
model for the evolution of substructure has been shown to 
be in reasonable agreement wi t h jV-body simulation results 
dTormen. Diaferio & Sver|[l998l: iTavlor & Babul| [200Tl [20031: 
iTaffoni et al.ll2003HZentner & Bullockll2003l) . 

When a halo merges and becomes a satellite, it is assigned 
an initial orbit in its host halo. The evolution of the orbit and 
the halo structure is then modeled including the effects of dy- 
namical friction, adiabatic orbital variations as the host halo 
grows, mass loss from the tides of the host halo at the orbital 
pericenter, and random orbital variations due to the gravita- 
tional deflections from other satellites. Satellites are also as- 
signed a core radius; when tides have disrupted the satellite 
down to the core radius, the satellite is considered to be com- 
pletely dissolved. For the application we are interested in this 
paper, the core radius of satellites is set to the black hole zone 
of influence whenever a black hole is present at the center of 
a satellite halo (the way halos are populated with black holes 
will be specified in §[5]): R z = GM./2ctq, where ao is the ve- 
locity dispersion of the Jaffe sphere at r = in the absence of 
a black hole and M, is the mass of a black hole at the cen- 
ter. We consider a satellite to be completely destroyed and 
merged with its host halo when Mj, < M., and we assume that 
their central black holes also merge at this point. If the satel- 
lite does not contain a black hole, its core radius is set to 1 % 
of the Jaffe radius. We describe these dynamical treatments in 
detail in the following subsections. 

3.1. Orbital Initial Conditions 

After a merger, a halo will initially move on an orbit with 
typical radius R m i, the radius of collapse of the host halo 
of mass Mh at the redshift of the merger introduced in §2.3. 
Therefore, we place the halo on an initial orbit with the same 
energy as a circular orbit at radius R co \: 



v 2 GM h 



In 1 + 



Rcoi 



(12) 



where we have used the potential of a Jaffe sphere with half- 
mass radius R^, and 



v 2 c (R co i) = 



GM h 
Rh+Rcoi 



(13) 



first generate a random radius r after having calculated the 
initial orbital energy eo with the distribution of equation (TBI) , 
and then we compute the potential energy at radius r and the 
modulus of the velocity vector at this radius required to have 
the orbital energy eo. The velocity is then assigned a random 
direction, which yields the orbital angular momentum, and the 
corresponding pericenter and apocenter of the orbit. 

3.2. Dynamical Friction 

We compute the rate at which the orbital radius of a satellite 
decreases owing to dynamical friction using the usual Chan- 
drasekhar formula. A satellite of mass M s at radius r s moving 
at speed v s in a field of particles with mass density ph(r s ) that 
move with a Maxwellian distribution of velocities of disper- 
sion a is subject to a dynamical friction acceleration equal to 



l6^lnAG 2 M sPh (r s ) 



^erf(X)-- e -* 2 
4 ' 2 



, (15) 



where X = v s /\/2a. Note that the exact value of adf in a 
Jaffe sphere is different because the velocity distribution is 
non-Gaussian and is described instead b y Dawson's integrals 
d Jaff d Il983t iBinnev & Tremaine 1987), but we neglect this 
small difference for simplicity. The Coulomb logarithm term 
is A = b mdx /b m i n , where Z? max is taken as the radius r„ and b m \ n 
is the larger of the satellite radius R s and GM s /v 2 (roughly 
the impact parameter at which the deflection angle would be 
1 radian for a point mass). 

This dynamical friction acceleration results in a loss of or- 
bital energy de/dt = — t^adf- Assuming the satellite moves on 
a circular orbit r c with v c (r c ) given eo, and using de/dr c = 
v c dv c /dr c + v 2 /r c , we find that the time derivative of the or- 
bital radius, f c , is 



d\nr c 



(16) 



which yields, after using equation (113t . 



f c _ -\6^\nkM s R h (r c + R h ) 
r c ~ P mt M h r c (r c + 2R h ) 



(17) 



In addition to the orbital energy, each halo is also assigned 
an initial orbital angular momentum. For this, we assume an 
isotropic velocity distribution, i.e., we assume that the phase- 
space density depends only on energy. This implies that the 
radial distribution of particles in orbits of energy eo, in the 
potential of the Jaffe sphere we have adopted, is 



P(r)droc4Trr dr / d v8(e— eo) oc r dr 




r+R h 



; +#h 
(14) 



where r max is the maximum radius at which the satellite can 
be located for the given initial energy eo. To choose the ini- 
tial orbital angular momentum of a halo after a merger, we 



where f or b = 2irr c /v c is the orbital period of the satellite. We 
use this equation in conjunction with formula (0 for the ve- 
locity dispersion at radius r c to evaluate the rate at which the 
orbit of every satellite decreases due to dynamical friction. 
For a more accurate calculation, this rate should be computed 
as a function of r s and the orbital eccentricity, but here we 
apply a constant rate of orbital energy decrease, independent 
of the eccentricity. We also require that the relative rate of 
change of the angular momentum per unit mass, L = r c v c , is 
independent of eccentricity: 



L _ rc Ue _ [c_ r c + 2flh 
L r c v c r c 2(r c +R h ) 



(18) 
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3.3. Adiabatic Orbital Variations 

As a halo increases its mass accreting new matter, the or- 
bit of a satellite will change in response to the time variation 
of the potential. We use adiabatic invariance to calculate the 
orbital changes due to the slow, gradual increase in the mass 
interior to the satellite orbit. Conservation of the action vari- 
able implies that the product r c v c must be conserved. Using 
the expression (fl~3T > for v c (r c ), we find that the orbital radius 
must vary at a rate 



on the fact that, to conserve the total energy of the system, 
the energy lost by satellites owing to the dynamical friction 
process must equal the energy gained by all the matter in the 
host halo (both smooth matter and satellites) owing to random 
encounters. We first compute the contribution to orbital per- 
turbations due to the dynamical friction of each z'-th subhalo. 
The orbital energy lost over a short time interval Af is 

AEi = EiAt = F dyn v c At . (22) 



r c + 2R h \ f c _ R h M h 
r c + R h J r c r c +R b M h 



(19) 



where Mh and Rh are the time derivatives of the mass and ra- 
dius of the halo. We also require that the angular momentum 
of the orbit is conserved to determine the eccentricity varia- 
tion. 

3.4. Tidal Mass Loss 

Subhalos are subject to mass los s due to tidal gravitational 
forces from t heir host halo (e.g., Bin nev & Tremain"elll987l: 
lGnedinll2003l) . It is useful to define the tidal radius r t , the ra- 
dius of the subhalo beyond which the tidal force from the host 
halo is stronger than the internal gravity of the subhalo. Par- 
ticles outside r t typically become unbound, resulting in mass 
loss. 

To model the tidal mass loss, we follow the simple prescrip- 
tion that each time a subhalo passes by its pericenter it loses 
all the mass that is outside the tidal radius r t . This tidal radius 
is obtained by requiring that the average density within r t in 
the subhalo and the average density of the host halo within the 
subhalo orbit are roughly the same: 



n = 



aM h (< D) 



1/3 



D 



(20) 



where D is the distance to the subhalo from the host halo cen- 
ter, and the dimensionless parameter is a = 2 for the Roche 
limit and a = 3 for the Jacobi limit. Hav ashi et al.l d2003l) 
investigated the tidal mass loss of dark matter subhalos in a 
static potential using iV-body simulations and found that the 
Jacobi limit is a good approximation. Here we adopt the Ja- 
cobi limit (a = 3). 

The mass and radius of the subhalo after each passage by 
its pericenter are modified according to 



M h =M h «r t )=M h 



n 



n+Rh 



R'h = 



n+Rh 



(21) 



This ensures that the velocity dispersion of the subhalo at 
small radius remains unaffected. 

3.5. Random Orbital Deflections 

Subhalos have their orbits perturbed by the gravitational in- 
teractions among themselves. These random perturbations 
may reduce the orbital pericenter, thereby hastening the de- 
struction of the satellite, or they may also increase the pericen- 
ter and allow the satellite to survive. Here we adopt a simple 
prescription to account for these orbital perturbations, based 



We note that the contribution from smoothly distributed mass 
of the host halo is negligible compared to the contribution 
from subhalos. Then, using the idea that the energy AEi is 
transferred to all the mass, AM,, of the host halo between 
the radii in which the orbit of the z'-th subhalo is comprised, 
the orbital perturbation Avj applied to the j-th subhalo is ob- 
tained by summing over all the contributions from the other 
subhalos, 




4 " 2 



,(23) 



where AM/ = Mh(< ^,,-)-Mh(< rj>;) is the mass of the host 
halo between the radii r^,, and r^,-, which are the apocenter 
and pericenter radii of the orbit of the z'-th subhalo. The sum- 
mation is made only over the subhalos whose orbits overlap 
with the orbit of j-th subhalo, and we adopt a simple weight- 
ing factor, 



Wl. 



minfaj, r AJ ) - max(z>, ; , r PJ ) 
r Aj-rpj 



(24) 



to account for the fraction of the time that the j-th halo spends 
in the region where it can interact with the z'-th halo. Although 
this time fraction is not linear with the radial width of the 
overlap region, we use this form for simplicity. Computing 
exactly the fraction of the energy that each satellite deposits 
in each radial shell would not introduce large differences on 
the statistical results in any case. 

In summary, at each interval Af , we add an orbital pertur- 
bation with amplitude Avj and arbitrary direction to the ve- 
locity of each subhalo, deflecting its orbit. We note that this 
method of perturbing the orbits of the satellites yields results 
that are independent of the timestep Af in the limit of a small 
timestep, as they should be. Because random perturbations 
in the velocity are added quadratically, the amplitude of the 
velocity perturbation at each timestep should indeed be pro- 
portional to VAf, as obtained in equation (|23| ). 

3.6. Transfer of Subhalos 

When a subhalo is being tidally disrupted, smaller satel- 
lites within the subhalo may be moved out of the system by 
the tidal forces, and have their orbit transferred to the larger 
halo within which the subhalo is orbiting. To incorporate this 
effect in our simple model, we remove satellites from their 
parent subhalos whenever their orbital radius is greater than 
the tidal radius r t , and we assign to them a new orbit in the 
larger halo. Since escapees follow the motion of their orig- 
inal subhalo, we assume that the new satellite orbit is close 
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halo mass (h~ l M Q ) 

FIG. 1 . — Halo mass function at z = 2. Solid curve: Halo mass function at 
Z — 2 in a high-density region that forms a halo of mass M h = IO^/t'Mq at 
z = 0. Dashed curve: Globally averaged halo mass function at z = 2. 

to that of their parent in the larger halo, although perturbed 
randomly by a velocity perturbation Av, which is applied at 
the position of the parent in the larger halo. The velocity per- 
turbation given to an escaping "sub-subhalo" is obtained by 
computing the velocity dispersion cr(r) at its position in its 
original subhalo and generating three one-dimensional veloc- 
ity perturbations along each axis from a Gaussian distribution 
with this dispersion. 

It happens occasionally that one of the satellites of a sub- 
halo acquires the escape velocity v^ir) = (2GMh/Rh)\n(l + 
Rh/r) owing to the random orbital perturbations. Then, it 
simply escapes the subhalo and is placed in a new orbit 
within the larger halo which is the same as that of its par- 
ent halo; in this case we add a random velocity perturba- 
tion with magnitude Av computed from energy conservation: 
Av 2 = v 2 (r)- v 2 sc (r). When this happens to a satellite of a 
parent halo that does not belong to any larger halo, the escap- 
ing subhalo is simply removed (in practice we find this does 
not occur very frequently). 

4. NUMERICAL MODEL 
4.1. Merger Tree 

We use the extended Press-Schechter formalism described 
in § 12. 1 I to generate the mass accretion histories of halos. We 
adopt a flat ACDM cosmology where the present Hubble con- 
stant is h = H /(100 kms^Mpc" 1 ) = 0.7, the matter density 
is Cl m = 0.3, the baryon density is fi/, = 0.02/T 2 , the power 
spectrum normalization is determined b y as = 0.9, and the 
primordial spectral index is n = 1 (e.g., ISpergel et al.ll2003t 
iTe^mark et al.l |2004). The po wer spectrum shape is o btained 
using the transfer function of lEisenstein & Hul (I1999I) . 

We are particularly interested in this paper in studying the 
origin of black holes (and the contribution of mergers to their 
growth) in massive galaxy clusters, such as the Coma clus- 
ter. We therefore generate a merger tree starting at z = with 
a halo mass Mh = \Q K~ 1 Mq, and identify all the progenitor 
halos with mass larger than the mass resolution M res . In order 



to accurately follow the growth and accretion history of halos 
with masses near M res , we actually generate the merger tree 
down to a smaller mass, which for a halo of mass M/,i is set to 
the minimum of 10~ 3 M/,i and M les . In this way, the variability 
of the halo mass growth due to the stochastic accretion of in- 
dividual halos is adequately reproduced. However, only halos 
with mass above M res are included when we follow their dy- 
namical evolution and the history of their central black hole. 

We adopt M res = 10 10 /z _1 M Q for the halos for which we fol- 
low the dynamical evolution. As we shall see in §5.1, the 
halos in which we place black holes in our model are of larger 
mass (^ 10 11 Ii~ 1 Mq) than our resolution limit, and the ha- 
los with no black holes are included in the calculation only 
for the purpose of taking into account the orbital perturba- 
tions they cause on the larger subhalos. We use a timestep for 
the merger tree corresponding to Az = 10~ 4 , which is small 
enough to ensure that P <C 1 in equation (0J. We have veri- 
fied the validity of our Monte-Carlo procedure by generating 
many merger trees for a large number of halos with the z = 
Press-Schechter mass function, and comparing the obtained 
mass functions at high redshift to the analytic Press-Schechter 
mass function. 

The halo mass function at z = 2 for a high-density region 
that is constrained to assemble into a halo of mass Mh = 
1O 15 /i~'M0 at z = is shown in Figure[T](solid line), compared 
to the globally averaged halo mass function at z = 2 (dashed 
line). Hereafter, we use the term "globally averaged" to refer 
to a halo mass functions averaged over the whole universe, 
without any restrictions to regions with special properties. As 
expected, massive halos are more abundant than average in a 
region selected for forming a massive cluster at a later epoch. 
The average of ten different realizations of a merger tree was 
used for the statistical results presented below in Figures [3] 
and |7] 

4.2. Examples of Dynamical Evolution 

It is useful to illustrate with a few examples the way that the 
dynamical processes described in §[3] work to determine the 
evolution of the subhalos, and ultimately the merger rates of 
black holes in our model. Figure[2]plots the redshift-evolution 
of the pericenter (r p ), apocenter (r^), Jaffe radius (R s ), and 
mass of six subhalos with black holes, which we have selected 
as illustrative examples (not randomly) from the merger tree 
for the formation of a halo with mass M h =10 15 /i- I M Q atz = 0. 
We generally refer to the largest halo present in the merger 
tree as the "cluster halo." 

In our model, a halo grows in mass from accretion and 
mergers as long as it merges only with objects smaller than 
itself, and it becomes a subhalo when it merges into a larger 
object. From that point on, the mass of the subhalo can only 
decrease as it experiences tidal disruption. In Figures and 
[2}?, the masses and radii of the subhalos are shown starting 
at the time they merge into the cluster halo. The decrease in 
mass and radius occurs at each pericenter passage, according 
to equation (f2~Tb . The thick solid and long dashed lines in Fig- 
ure |2jz show the mass of the cluster halo and of the central 
black hole, while the radius of the cluster halo is shown as the 
thick solid line in the other three panels. 

The evolution of the orbits is shown in Figures |2j: and|2j/, 
which plot the pericenter and apocenter as a function of time. 
There is an average tendency for the orbits to shrink because 
of dynamical friction. The rate at which the orbits shrink in- 
creases with subhalo mass and with proximity to the center. 
As the subhalos approach the center, their masses and radii 
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FIG. 2. — Examples for the dynamical evolution of six subhalos hosting a black hole, after they merge into the large cluster halo that reaches a mass 
M = 1O 15 /t'M0 at z = 0. The thin lines in the four panels show the subhalo mass, radius (Rs), pericenter (rp) and apocenter r^) of the six satellites as a function 
of redshift. In the upper left panel, the thick curves are the parent halo (solid) and black hole (dashed) mass. The thick solid curve in the other three panels is the 
luster radius (i?h). 



are reduced by tidal disruption. It should be noted here that 
when the mass of both the subhalo and cluster halo inside ra- 
dius r vary linearly with r (as is the case for Jaffe spheres in 
the inner parts), the mass of a subhalo after tidal disruption 
at pericenter r p will be proportional to r p , and the ratio of 
the subhalo mass to the cluster halo mass enclosed within the 
subhalo orbit remains constant. Therefore, the ratio of the dy- 
namical friction time to the orbital time remains constant also, 
and the subhalo continues to spiral toward the center in a con- 
stant number of orbits in each logarithmic radial interval. The 
subhalo is therefore completely destroyed in a finite time, and 
the black holes are then assumed to merge. In the example in 
the figures, two halos are destroyed in the center at z — 0.56 
and z — 0.25, causing the jumps in the mass of the central 
black hole. The mass ratios of the two mergers are 30% and 
83%, respectively. There is one more black hole that merges 
at z — 0. 16 with a mass ratio of 168% (not shown in the Figure 
for the sake of clarity). Another four halos are shown which 
are gradually decreasing their orbital radius but still survive 
by z = 0. The halos added to the cluster at z — 0.25 - 0.3 have 
a very slow rate of orbital decay, and their pericenter is large 
enough to avoid any tidal disruption. Their apocenters and 
pericenters show random variations owing to the gravitational 



encounters among subhalos. 

In general, subhalos are destroyed when the dynamical fric- 
tion time is shorter than the age of the system. Subhalos that 
avoid tidal disruption have an initial dynamical friction time 
long compared to the age of the universe. If the cluster halo 
were to remain perfectly static, eventually all subhalos would 
spiral to the center given a sufficiently long time. However, 
this does not happen because as the cluster halo continues to 
grow in mass, large subhalos that are continuously merging 
cause random perturbations on the orbits of smaller subha- 
los which can increase their orbital radii and therefore their 
dynamical friction times. Additionally, whenever the halo 
within which a satellite is orbiting merges into a larger sys- 
tem and is tidally disrupted, the satellite is placed in an orbit 
in the larger system with a much longer dynamical friction 
time. In this way, small subhalos can survive indefinitely or- 
biting in larger halos as long as the halos continue to merge 
and grow. 

The important dynamical effects determining the rate of 
black hole mergers obtained in our model are the fraction of 
satellites that can spiral all the way into the center of their 
parent halo by dynamical friction before they are scattered, 
and the rate at which they are tidally disrupted as they pass 
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FIG. 3. — Solid line: Number of subhalos found within the IO^/t'Mq 
halo at z = with a mass fraction above a value M s /M/,, obtained in our 
model from the average of 10 numerical realizations. Shaded region: sta- 
tistical uncertainty, computed from the error on the mean of the ten realiza- 
tions. Dashed line: Analytic subhalo mass f unction found to fit resu lts of full 
numerical simulations of halo substructure I Vale & Ostriker 2006). Bottom 
panel: Percentage deviation of our calculation relative to this fit. 

near the halo center. Hence, the black hole merger rates we 
present in §6 are most sensitive to the rate of dynamical fric- 
tion and mass loss by tidal disruption as a function of orbital 
eccentricity and semimajor axis, as well as the rate of orbital 
perturbations by other satellites. 

As a test of our analytic model, we compare the subhalo 
mass function we obtain in the cluster halo at z = 0, plotted in 
Figure[3] with results that have been obtained from numerical 
simulations. We include only subhalos on orbits with a cir- 
cular radius smaller than the cluster halo radius, to compare 
with numerical simulation results where subhalos are identi- 
fied as bound objects inside the virial radius of the parent halo. 
Our subhalo mass function is computed from the average of 
ten realizations of the merger tree of the cluster halo. We fit 
the result to a Schechter mass function for the subhalo mass 
function, 

/ M s \~ a ( M s \ dM s 
N(M s \M h)dMs =A^ w J exp^-_J_, (25) 

where a is the power- law slope and /3Mh is the cutoff mass 
(IVale & Os triker 2006). The normalization is given by 



r[2-a]-r[2-a,M s , max //3M h ] ' 

where the total mass fraction in subhalos is 7, M s . max = 0.5Mh 
is the maximum subhalo mass, and T[x] and r[a,jc] are the 
Gamma and the incomplete Gamma functi ons, respectively . 
We fix the cutoff mass fraction to [3 = 0.3 dShaw et alj|2005l) 
and fit the analytic function to our results, allowing the param- 
eters a and 7 to vary. The resulting best-fit subhalo mass func- 
tion is shown as a dashed line, and has values a = 1.82 and 
7 = 0.12, in good agreement with typical valu es found in nu- 
merical simulations, a = 1 .7 - 1 .9 and 7 ~ 0. 1 dDe Lucia et alJ 
I2004t IShaw et alJl200l IZentner et al.ll2005h . Our model ba- 
sically agrees with the results found in simulations within the 
uncertainty. 



5 . BLACK HOLE MASS FUNCTION IN CLUSTERS AT Z = 2 

This section describes the way we initially populate halos 
with central black holes. Our dynamical evolution model de- 
scribed previously tells us the rate at which halos merge all the 
way to their central cusps, and we assume that the merger of 
the central black holes follows immediately after the merger. 
The simple model that we now describe for the initial popula- 
tion of black holes in halos allows us to derive a rate of black 
hole mergers from the rate of halo mergers that proceed all 
the way to their centers. 

A full model of the evolution of the black hole mass func- 
tion including the effects of gas accretion and mergers would 
include the gradual growth of black holes from gas accre- 
tion as constrained by observations of the quasar luminosity 
function at each redshift, and would simultaneously follow 
black hole mergers from the dynamical evolution of host ha- 
los. Here we make a simplifying assumption that separates the 
two evolutionary factors for the growth of black holes: gas ac- 
cretion and mergers. Observationally, we know that most of 
the radiative energy of luminous quasars was emitted over a 
narrow range o f reds hift, 1 .5 < z < 4 (e.g JCroom et~aT]l2004t 
iRichards et alJ l2006h . Hence, we make the approximation 
that gas accretion creates all the black holes instantaneously 
at z = 2, roughly at the epoch of the peak in quasar activity, 
and that subsequently the black hole mass function evolves 
only by mergers. We choose to place the black holes in halos 
at z = 2 with their observed mass function at z = 0. If black 
hole mergers cause relatively small changes in the black hole 
mass distribution, then the evolved mass function should still 
be approximately correct. 

The black hole mass function at z = is often estimated 
by using observed luminosity functions or velocity disper- 
sion measurements of galaxies, and observed correlations 
with the central black hole masse s. Here we use the black 
hole mass function obtained by ISteed & W einberg ( 120031) . 
This is derived fro m the distribution o f early-type galaxy ve- 
locity dispersions dSheth et al. 2003) using the relation be- 
tw een the black hole m ass M. and velocity dispersion a 
in iTremaine et alJ ([2002), and adding also a contribution 
from spiral galaxies (these become the dominant host galax- 
ies of black holes with M , < 10 8 M Q , using the estimate of 
lAller & Richstonel (120021) ). We also include an intrinsic scat- 
ter to the M. - a power-law relation of 0.5 dex, with a log- 
normal distribution in black hole mass. 

We assign black holes to host halos at z = 2 by assuming 
that all halos that have not merged into any larger object con- 
tain a black hole by z = 2, and that there is a monotonically 
increasing, one-to-one relation between the halo mass, M n , 
and the corresponding black hole mass, M.(Mh), at this red- 
shift. This implies that the cumulative number densities of 
halos with mass above Mh must be equal to the cumulative 
number density of black holes with mass above M.(Mh). This 
is illustrated in Figure Hh, where the globally averaged halo 
and black hole mass functions are plotted. As an example, 
halos of mass M h = IO'-V'Mq have the same number den- 
sity as black holes of mass 6 x 10 8 M Q , so we place a black 
hole of this mass in each halo of mass Mh. We apply this re- 
lation down to a minimum black hole mass of 10 6 M Q , which 
corresponds to a halo mass Mh > 8 x lQ ll h~ l M Q . 

Note that the merger tree and the dynamical evolution of 
halos are computed from z ^> 2, even though black holes are 
only assigned to halos at z = 2. Furthermore, the black holes 
have no impact on the dynamical evolution of their host halos 
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FIG. 4. — Populating halos with black holes at z = 2. 
functions at z = 2. As an example, black holes of mass • 
(b) Solid curve: black hole mass function in the cluster region at z = 2; dashed curve: global black hole mass function at z = 2. 



(a) Assumed relation of the halo and black hole masses, determined from their cumulative global mass 
6 X 1O 8 M0 are placed in halos of mass ~ IO^/t'Mq by matching their cumulative number densities. 



in our model. When orbiting halos merge all the way to their 
centers, we simply add the mass of their two black holes to 
obtain the mass of the new central black hole. 

In Figure @}?, the black hole mass function in the cluster 
halo at z = 2 obtained in this way is compared to the glob- 
ally averaged black hole mass function. The black hole mass 
function in the cluster halo that we simulate is of course not 
the same as the globally averaged black hole mass function 
because, owing to the biased spatial distribution of halos, the 
region that collapses into a massive halo at z = contains an 
excess density of halos at z = 2 and therefore an excess den- 
sity of black holes compared to the global average. The mass 
function is shown as the number of black holes per unit of co- 
moving original volume before the cluster collapses. There- 
fore the larger number density of black holes is not due to 
the physical collapse of the cluster halo, and reflects only the 
biased distribution of black holes. 

6. IMPACT OF MERGERS ON BLACK HOLE GROWTH 

In this section, we present the results on the impact of merg- 
ers on the black hole mass function in a massive cluster and 
the way black holes are distributed in a cluster. We do not 
consider any possible ejections of black holes due to gravi- 
tational wave recoil, and we neglect the mass lost due to the 
emission of gravitational waves during mergers. 

First, we examine the results of one particular realization of 
a 10 15 /i -i Mq halo. The left panel in Figure|5]shows the initial 
and final values of the black hole mass and the host subhalo 
mass for all the black holes hosted by satellite halos of the 
cluster. The small dots are the values at z = 2, and they follow 
the black hole mass-halo mass relation that we impose as ini- 
tial conditions. At z = 0, the black hole and halo masses are 
shown as filled symbols. The central black hole in the cluster 
is the one hosted by the main halo with mass 1O 15 /i _1 M0, and 
all other halos have become satellites in the cluster. Subhalos 
with their own satellite systems that have survived tidal dis- 
ruption by the main cluster halo are shown as filled squares, 
and satellites of subhalos are shown as triangles. Note that 
the evolution of the halos from z = 2 to z = may either in- 
crease the halo mass due to mergers and accretion of other 



halos, or decrease the halo mass due to tidal stripping and 
disruption when they become subhalos in a larger halo. On 
the other hand, the black hole mass can only increase if two 
black holes merge, following a complete disruption of the host 
subhalo. To facilitate the identification of progenitor black 
holes, we connect filled symbols (z = 0) and dots (z = 2) of 
the black holes that have experienced any mergers. The right 
panel in Figure [5] shows the final orbital radius of each satel- 
lite in the cluster and the black hole mass, showing little cor- 
relation between these two quantities, except perhaps at the 
highest masses where the black holes may tend to be hosted 
by satellites closer to the cluster center. 

Most of the black holes with masses above - 10 75 M Q ex- 
perience mergers, and some of them increase their mass in a 
substantial way. The central black hole is the one that grows 
the most, experiencing many mergers and growing its mass by 
a factor ~ 5. At lower masses, most black holes do not expe- 
rience any merger; for these, the filled symbol is at the same 
black hole mass as the small dot (note that some small dots do 
not have any corresponding large symbol because they have 
merged with a more massive black hole). In our dynamical 
model, the satellites representing the host galaxies of these 
black holes form galaxy groups and later become part of a 
larger cluster, where their dynamical friction time is too long 
to produce any significant decay. Although the outer parts of 
these satellites may be tidally disrupted (producing the large 
reductions in halo mass for most of the cases shown in the 
figure), the core of the satellite containing the black hole sur- 
vives and remains as a satellite halo. We need to caution that 
the rarity of mergers for low-mass black holes may be caused 
in our model by the fact that we do not include black holes 
with initial mass below 10 6 M©, and we place all the black 
holes at z = 2 in field halos (halos that are not satellites in a 
larger halo). 

The black hole that ends up in the cluster center is identified 
as the one that was present in the original halo which, merging 
only with other halos smaller than itself, became the cluster 
halo at z = 0. In any merger-tree simulation of the formation of 
a cluster, there is one and only one initial halo that has always 
merged with other halos smaller than itself, which is found 
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FIG. 5. — Distribution of satellite halos containing a nuclear black hole in one realization of the cluster halo. Left panel: Filled symbols indicate masses of 
black holes and their host subhalos at z = 0. Small dots give the initial relation between halo mass and black hole mass at z = 2. For black holes that have merged, 
the initial and final masses are connected; for other black holes, only the halo masses have changed. The asterisk and open circle identify the central halo and 
black hole in the cluster. Squares indicate subhalos that have their own satellite halos with black holes, and these satellites are shown as triangles. Right panel: 
Orbital radius of satellites hosting a black hole in the cluster (the central black hole is not plotted). The horizontal bar is the radius of the cluster. 

elliptical galaxy (or the bulge component of a spiral galaxy). 
This observed relation tells us that, irrespective of whether a 
giant galaxy is considered to be the central one in a cluster or 
a satellite, its stellar mass should follow the same correlation 
with black hole mass that was imprinted during the quasar 
epoch at high redshift. Mergers will increase the mass of a 
black hole and the spheroidal stellar component by the same 
factor (neglecting gravitational wave emission and ejected 
stars), so in any clusters where the central black hole is not 
the most massive one, it should probably also be true that the 
central galaxy is not the one with the largest stellar mass, and 
that the most massive black hole is in the galaxy with the most 
massive stellar component. 

Figure [7] summarizes the results of the ten realizations. The 
initial halo mass functions at z = 2 of the ten realizations are 
shown in the first panel, together with the average mass func- 
tion. The mass functions in the realizations contain many 
more massive halos than average because of the bias intro- 
duced when selecting a region that forms a IO^/t'Mq halo 
at z = 0. The second panel shows the black hole mass func- 
tions. For the same reason, the initial mass function contains 
many more massive black holes per unit volume than the aver- 
age mass function. Black hole mergers result in a substantial 
change in the black hole mass function, in particular a large 
increase in the number of black holes of M > 5 x 10 9 M Q . 

Because of this increase, our initial conditions in which the 
z = 2 black hole mass function matches the z = mass function 
are not fully consistent with our results. However, one should 
note that Figure [7j? shows the change in the black hole mass 
function in a 10 15 /! _1 M Q cluster, and that such clusters are 
rare; the global evolution of the high end of the mass function 
will be weaker. 

In the last two panels, thick solid histograms show the num- 
ber of mergers experienced by each black hole in terms of the 
black hole mass, and the average factor by which the black 
hole mass increases as a result of these black hole mergers. 
Figure 6c shows that most black holes of M > 1O 8 M expe- 
rience mergers. However, their mass growth is small except 
for the most extremely massive black holes, M > 3 x 10 9 M Q . 
The majority of black holes, with lower masses, grow only 
by a small factor because they often merge with black holes 
of much lower mass than themselves. At low masses, black 



by tracing back the merger history from the final cluster and 
choosing always the most massive progenitor halo at every 
merger event. This halo is not necessarily the most massive 
halo at z = 2. Because the growth rates of halos are stochastic, 
there may be halos of high mass at z = 2 which grow slowly, 
while other halos of originally lower mass grow faster and 
become more massive by the time that all the halos merge to 
form the final cluster at z = 0. This is in fact the case in the 
example of Figure [5] where the initial halo that becomes the 
center of the final cluster is identified by an asterisk. This 
central progenitor has a mass of ~ 10 12 7 /i _1 Mq at z = 2 and it 
grows rapidly to become the cluster halo, while several other 
halos that were originally more massive (the largest having 
nearly 10 i4 !i~ i Mq at z = 2) become satellites as they merge 
at late times with the fast-growing halo. When this situation 
occurs, the black hole in the central galaxy will not be the 
most massive one. In the present example, the black hole in 
the cluster center would have a mass of only 10 9,2 M Q , while 
black holes in other satellite galaxies have masses as high as 
~ 10 98 M Q . 

This situation, however, is not the most common one, as 
can be seen in Figure|6l where an additional nine random re- 
alizations of the formation of a lO 15 /i _l M halo are examined. 
In the majority of cases, the halo that becomes the final clus- 
ter at z = was the most massive halo at z = 2. Out of the 
ten random realizations we have performed, there is only one 
case where black holes in satellite galaxies are substantially 
more massive than the black hole ending up at the center, but 
there are several where black holes in satellite galaxies are of 
comparable mass to the central one. We also see that there 
is a substantial dispersion for the mass of the central black 
hole in a cluster of fixed mass (10 15 /z -1 Mq), going from val- 
ues as high as 1.5 x 10 10 M Q to as low as 2 x 1O 9 M . If we 
take the most massive black hole in the cluster (rather than the 
central one), then the low-end of this mass range increases to 
4 x 10 9 M Q , still implying a substantial dispersion. 

In practice, the distinction between the black hole in the 
cluster center and black holes in satellites is ambiguous ob- 
servationally. Clusters typically have substructure and may 
have more than one galaxy that could qualify as "central." 
However, observationally we know that there is a linear cor- 
relation between black hole mass and the stellar mass of an 




FIG. 6. — Distribution of the black holes in the nine additional realizations of the merger tree of the cluster halo. The symbols are the same as Fig. [5] 



holes are removed by merging into more massive objects, 
causing a decrease in the mass function of a factor of ~ 2 
(FigJTJ?). We are cautious in interpreting this result because 
it may be severely affected by our approximation of placing 
the black holes initially at z = 2: had we followed the evolu- 
tion of black holes from a higher redshift, when many more 
low-mass halos are present, many more mergers of low-mass 
halos (and their central low-mass black holes) would proba- 
bly have occurred. The minimum black hole mass of 1O 6 M 
in our simulation should probably also be increased in order 
to correctly follow the black hole mergers at higher redshift, 
and a model of black hole growth from gas accretion that is 
distributed over redshift (instead of occurring instantaneously 
at z = 2, as we have assumed here) would need to be included. 



A fundamental process affecting the way black hole merg- 
ers can occur in our model is the mechanism of random orbital 
perturbations described in § 3.5 to take into account the in- 
teractions of the satellite halo containing the black hole with 
other satellites. If these perturbations are absent, a satellite 
can only continue to undergo orbital decay by dynamical fric- 
tion, unless its parent halo merges into another object. The 
perturbations allow a satellite to randomly gain orbital energy 
and in this way avoid tidal destruction. It is therefore useful to 
check how much our merger rates are altered if we suppress 
these interactions, to see if the results are highly sensitive to 
their presence and to the detailed way in which we implement 
the perturbations. This is quantified in Figure 7, in panels 
(c) and (d), where thin solid histograms show the results with 
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FIG. 7. — (a) Thin solid lines: Halo mass function in the M = 10 A Afg cluster at z = in each realization. Thick solid: Halo mass function averaged over 
the ten realizations. Thick dashed: Globally averaged halo mass function, (b) Black hole mass function in the cluster initially (at z = 2) and at z = 0, after the 
mergers have occurred. The thick dashed line is the globally averaged black hole mass function, (c) and (d) Individual small dots are the number of mergers and 
mass growth factors for black holes that have merged in ten realizations. Thick solid lines are averages. Thin solid lines are the number of mergers and mass 
ratios, when interaction between subhalos is ignored. 



the perturbations turned off. As expected, the number of black 
hole mergers increases, although by a relatively small amount, 
while the mass growth remains similar. Hence the results are 
not greatly affected by the presence of the perturbations. We 
also find that the enhanced disruption of subhalos when per- 
turbations are turned off results in a subhalo mass function for 
which only 7% of the total mass is contained in subhalos, less 
than is found in numerical simulations. Thus, our physically 
motivated recipe for including orbital perturbations produces 
better agreement between the analytic model and numerical 
simulations (see Figf3]l. 

7. CONCLUSIONS 

We have developed a dynamical model of substructure in 
dark matter halos based on numerical merger trees and an- 
alytic models of the dynamical evolution of satellites, with 
the purpose of modeling the rates of mergers of nuclear black 
holes. We have improved on previous work by tracing the 
evolution of satellites in a satellite halo, and by accounting 
for satellite-satellite interactions. Our model reproduces the 
analytic subhalo mass function found previously in numerical 
simulations, forM, > 10" 5 M h in a halo of M h = 10 15 /z _1 M Q . In 



this paper, we have focused on examining the effect of merg- 
ers on the abundances of the most massive black holes in clus- 
ters of galaxies, using a simplified model where black holes 
have completed their growth by gas accretion by z = 2 and re- 
main with a fixed mass after that, unless they undergo merg- 
ers when their host galaxies merge. We have assumed that the 
merger of two satellites in a halo always leads to the merger of 
their central black holes (which maximizes the merging rate 
for black holes). 

We find that mergers have an important impact on the most 
massive black holes in the cluster. For the most massive black 
holes in a cluster, with M > 10 9 M©, a growth in mass by a fac- 
tor ~ 2 is typical. Our ten realizations of a \0 15 IT 1 Mq cluster 
include four black holes with a final mass of 1 - 1 .5 x 10 10 M Q , 
but without mergers such black holes would be extremely 
rare. 

Mergers are less important for black holes of lower mass. 
The main effect found in our simulation on low-mass black 
holes is that they are depleted by a factor ~ 2 because of their 
mergers into more massive black holes. This result may be af- 
fected by our procedure of initially populating the halos with 
black holes at z = 2, and by the minimum black hole mass of 
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10 6 M© that we use. We also caution that the effect of de- 
pletion in the global black hole mass function may be much 
weaker than the effect in the cluster environment, because by 
the epoch of z = 2, most of the low-mass black holes that were 
formed in a region collapsing in a massive cluster at present 
have already merged in larger objects. 

We have also found that the most massive black hole is not 
always located in the central galaxy of the cluster, if we de- 
fine the central galaxy to be the one that formed in the halo 
that has always merged with other halos smaller than itself. 
In practice, it may be difficult to assign observationally which 
galaxy in the cluster corresponds to this central one; more- 
over, the linear relation between black hole mass and bulge 
stellar mass should not be altered. 

Our model may overestimate the importance of mergers be- 
cause we have assumed that, once the nuclei of two satel- 
lite halos have approached each other down to the radii of 
influence of the black holes, the black holes will always 
merge. Actually, two black holes may become stalled af- 
ter they have scattered all the stars interior to their orbit 
as well as any remaining stars in the loss-cone, and before 
reaching the radius of orbital deca y by gravitational waves 
(Begel man. Rees & Blandford 1 120061) . so it may be that not 
all black hole binaries merge. It is unlikely that a majority of 
mergers end up in stalled binary systems, however, because 
few binary black holes in galactic nuclei are found, and be- 



cause the presence of gas is likely to drive the binary to a 
merger in any case (Gould & Rix 2000). Moreover, models 
that try to estimate the probability of ejection of black holes 
during three-body interactions find that only a small fraction 
of the black holes can be lost in this way (Volonteri et al. 
2003). 

Therefore, the results of our paper support the idea that 
black hole mergers are important in shaping the black hole 
mass function at low z, particularly affecting in a strong way 
the abundance of the most massive black holes in rich clusters 
of galaxies. In a subsequent paper, we shall examine more 
generally the impact of mergers on the evolution of black 
holes outside of massive clusters. 
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